# This code replicates Figure 9
rm(list=ls())
setwd("")
library("gridExtra")
library(foreign)
library(ri)
library(ggplot2)
data <- read.dta("data.dta")
data <- data[complete.cases(data$exp1_leave), ]

#### Adventist
y <- data$exp1_leave
Z <- data$treatment_adventist
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE

## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

p <- dispdist(distout, ate)[4]
ate
distout <- as.data.frame(distout)
m1 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey") + 
      geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
      annotate("text", x = -0.18, y = 42, label = "ATE: -0.135 \ \ \ \ ") +
      annotate("text", x = -0.18, y = 40, label = "P-Value: 0.048") +
      labs(title = "Adventist Treatment Effect") + 
      labs(x = "Estimated ATE") +
      labs(y = "Frequency")
m1


#### Mixed

y <- data$exp1_leave
Z <- data$treatment_mixed
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE
## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

dispdist(distout, ate)[4]
ate
distout <- as.data.frame(distout)
m2 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey") + 
  geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
  annotate("text", x = -0.18, y = 42, label = "ATE: -0.003 \ \ \ \ ") +
  annotate("text", x = -0.18, y = 40, label = "P-Value: 0.359") +
  labs(title = "Mixed Treatment Effect") + 
  labs(x = "Estimated ATE") +
  labs(y = "Frequency")
m2



#### Peruana

y <- data$exp1_leave
Z <- data$treatment_peruana
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE

## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

dispdist(distout, ate)[3]
ate
distout <- as.data.frame(distout)
m3 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey") + 
  geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
  annotate("text", x = -0.18, y = 42, label = "ATE: 0.091 \ \ \ \ \ ") +
  annotate("text", x = -0.18, y = 40, label = "P-Value: 0.123") +
  labs(title = "Peruana Treatment Effect") + 
  labs(x = "Estimated ATE") +
  labs(y = "Frequency")
m3


#### Maranatha

y <- data$exp1_leave
Z <- data$treatment_maranatha
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE

## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

dispdist(distout, ate)[3]
ate
distout <- as.data.frame(distout)
m4 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey") + 
  geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
  annotate("text", x = -0.18, y = 120, label = "ATE: 0.133 \ \ \ \ \ \ ") +
  annotate("text", x = -0.18, y = 113, label = "P-Value: 0.022") +
  labs(title = "Maranatha Treatment Effect") + 
  labs(x = "Estimated ATE") +
  labs(y = "Frequency")
m4


grid.arrange(m1, m2, m3, m4, ncol=2, nrow=2, 
             widths=c(3, 3), heights=c(3, 3))
pdf("obedience_randomization.pdf", width = 9, height = 9)
grid.arrange(m1, m2, m3, m4, ncol=2, nrow=2, 
             widths=c(3, 3), heights=c(3, 3))
dev.off()

